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In this article, we develop a numerical method to find optimal control pulses that accounts for the 
separation of timescales between the variation of the input control fields and the applied Hamilto- 
nian. In traditional numerical optimization methods, these timescales are treated as being the same. 
While this approximation has had much success, in applications where the input controls arc filtered 
substantially or mixed with a fast carrier, the resulting optimized pulses have little relation to the 
applied physical fields. Our technique remains numerically efficient in that the dimension of our 
search space is only dependent on the variation of the input control fields, while our simulation of 
the quantum evolution is accurate on the timescale of the fast variation in the applied Hamiltonian. 
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I. INTRODUCTION 

Quantum information promises to significantly speed 
up hard computational tasks [1]. There is a range of 
candidate systems for quantum computing hardware in- 
cluding atomic systems [2], spins [3], and solid-state cir- 
cuits [4] . Quantum algorithms are typically implemented 
in these as a sequence of quantum logic gates, each of 
which is implemented as a short temporal pulse of an 
external field. These pulses have to be adapted to the 
properties of the hardware at hand, which is often im- 
perfect. This is the task of control theory — finding an 
ideal pulse shape that implements the desired gate (or 
state transfer). 

Control theory is rooted in applied mathematics and 
has been successfully mapped to quantum physics. The 
first applications of quantum control were in physical 
chemistry [5, 6] and nuclear magnetic resonance (NMR) 
[7, 8], but it has since been applied to many diverse sys- 
tems such as atoms [9], superconducting qubits [10, 11] 
and semiconductors [12] to name a few. There exists a 
wealth of analytical and numerical methods [13-16, 18- 
20], the latter most often employing the GRAPE [13] or 
Krotov [14, 15, 19] algorithms. 

Novel challenges arise in engineered solid-state quan- 
tum systems, such as superconducting qubits, where 
quantum coherence is typically limited to microseconds, 
necessitating short pulses [21-25]. Also, in these systems, 
the control design is limited to the microwave frequency 
range where current state-of-the-art electronics restrict 
the variation in controls to a few nanoseconds [24, 25]. 
Thus, the techniques that arc otherwise well-developed in 
other contexts need to be tailored to these constraints. 
Properly quantifying the shape and controllability of fast 
pulses using slow modulation is an important component 
of lowering effective error rates. There has been some 
progress in this direction with the development of opti- 
mal control methods such as DRAG [10, 26] and Rcfs. 
[11, 27]. An added benefit is that faster gates linearly 



speed up computation in these systems. 

Most numerical control methods assume that the sam- 
pling rate of the pulse shape that can be optimized is 
identical to the sampling rate of the control field; in other 
words, either the control fields can be fully shaped in 
real time, or we shape the envelope of a driving field 
represented by the Hamiltonian in an appropriate ro- 
tating frame where all remnants of time-dependence on 
the scale of the driving field disappear. In both cases, 
the time-evolution across a step of the shaped pulse, 
a pixel, can be straightforwardly approximated using a 
time-independent Hamiltonian in either frame. 

In this article, we are addressing a more general 
case in which the Hamiltonian across the pixel is time- 
dependent. This can occur for a wealth of reasons. If the 
piilse-shaping hardware has an internal filter with a band- 
width that is much smaller than the input control sample 
rate, then the control fields become smooth rather than 
a sequence of plateaux. Secondly, for many applications 
where one requires a high-fidelity operation, performing 
a rotating wave approximation (dropping fast rotating 
terms from a Hamiltonian) becomes invalid when the 
control pulses are extremely short. That is, the coun- 
tcrrotating terms should be taken into the integration of 
the Schrodinger equation even though the input controls 
can only be changed on a time scale which is much slower 
than these terms. More generally, if multiple Fourier 
components of a driving field arc applied in order to, e.g., 
induce AC Stark- and Zceman-shifts, not all of these fre- 
quencies can be eliminated by a suitable transformation. 

To account for this effect, we introduced an extra level 
of discretization. We separate the discretization neces- 
sary for integrating the quantum evolution from that of 
the discretization of the control parameters . We show 
with some examples that only a few extra subdivisions 
are needed to greatly increase the accuracy of the op- 
timization. Moreover, since the number of controls re- 
mains the same, the search space of possible pulses does 
not change and the time-cost of the algorithm is only 
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linearly affected with the number of subdivisions. 

The paper is organized as follows: in Sec. II, we out- 
line the control problem for optimizing quantum opera- 
tions, describe conventional approaches, and discuss their 
limitations; in Sec. IV, we establish our coarse-grained 
approach to address these limitations; specific transfer 
functions for piecewise cubic splines, filters, and Fourier 
components are given in Sec. V; this methodology is ap- 
plied to select problems in section VI and concluded in 
Sec. VII. 



II. THE CONTROL PROBLEM 

In this method the standard control problem is to find 
a way to implement a desired unitary operation L^idoab 
given discrete controls each of length N. It is standard 
to assume that the underlying Hamiltonian is of the form 



(2.1) 



where Hk is the control Hamiltonian and Ck{t) is a con- 
tinuous field which is parameterized by the controls u^,. 
The goal is to choose controls such that at the final 
time t = T , the dynamics of this system approximates 
the target unitary operator: J7ideai- To find these controls 
one first needs to simulate the dynamics of the system, 
either analytically or numerically, given the Hamiltonian. 
In addition, a cost function C is needed to describe how 
close a simulated operation is to J7ideai- This allows us to 
evaluate the performance of each and thus search the 
space of controls for an optimum. 

The most general measure for comparing two quantum 
operations would be the diamond norm [28], 



C= ||PQ(Mdcal-A^(T))PQ||o. 



(2.2) 



where (T) is our simulated quantum operation and Pq 
is a projector onto the subspace of interest. The ideal 
quantum operation is typically a unitary, formally writ- 
ten in Liouville space as ideal — CAdcai ® C^idcai where 
C^ideai IS a complcx conjugate. However, in the interest 
of defining a more tractable search problem, we gener- 
ally ignore the incoherent dynamics and instead require 
that the controls are time-optimal, i.e. shorter control 
sequences imply more coherence. There are exceptions 
to this rule for systems with non-Markovian decoher- 
ence [29] or approximately decoherence-free subspaces 
[15, 30, 31]. Minimizing the diamond norm is generally a 
hard problem; however, for many cases, minimizing the 
average error is sufficient. This can be done using the 
Frobenius-norm as the cost function, which for unitary 
maps is equivalent to maximizing [13, 32] 



* = ilTr(^/LalC^(^)PQ)l^ 



(2.3) 



The mathematical statement of the problem is to op- 
timize $ with respect to the vector u^. It is not im- 
mediately apparent how one should perform this opti- 
mization. One of the simplest approaches is the steepest 
ascent (or gradient search) algorithm. If we consider the 
multi-dimensional surface (control landscape) formed by 
^(ufc), steepest ascent is an iterative update procedure 
that locally examines the landscape at each iteration and 
provides a new by moving in the direction that in- 
creases $ the greatest. We take some initial guess for 
the controls and iteratively update according to the rule 



Ufe ^ Ufc + eVu'fc$(ufe) 



(2.4) 



where e is a unitless step matrix (such as the inverse 
Hessian for BFGS-type algorthims). Given an arbitrary 
initial configuration for the control fields, the algorithm 
follows a steepest ascent to a local optimum, at least 
in the case of small e. For simplicity, in this work we 
take e to be a scalar which is chosen adaptively. It is in 
no way clear that this procedure will produce anything 
other than local maxima, but remarkably, for the types 
of $'s we examine in quantum control, the landscape 
is sufficiently under-constrained that gradient searches 
find global maxima with high probability. This has deep 
implications for the topology of the control landscape 
[33, 34]. 



III. THE STANDARD ALGORITHM 

The GRAPE algorithm [13] is an example of a gradient 
search technique. In order to numerically integrate the 
Schrodinger equation it is necessary to discretize the evo- 
lution, normally into piecewise constant intervals. The 
approach traditionally used in the GRAPE algorithm is 
to also use this discretization as the map between the 
continuous fields and the control vector. That is 



N~l 

Cfe(i) = Uk,j V^J (i, Ai), 



(3.1) 



with the rectangle function 



(t, A<) = [Q[t - j^t) - Q[t - {j + l)At)] , (3.2) 

where O is the Heaviside unit step function. Evaluat- 
ing the Schrodinger evolution using these fields and the 
Hamiltonian in Eq. (2.1) yields 



U{0,T)^ H U„ 



where 



exp 



-iH{jAt)At 



(3.3) 



(3.4) 



where dn is the dimension of the computational subspace. 



the product running backwards is to enforce time- 
ordering, and NAt = T. 
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A principle insight of GRAPE is that the gradient of 
$ can be computed more efficiently than that of an ar- 
bitrary function. Specifically, the derivative of $ with 
respect to one of our control variables is 



duk,j 



duk,j 

X TV (C/idealPQC/(T)tPQ) 



(3.5) 



where 




(3.6) 



The speedup over gradient searches based on a naive ap- 
proach where one queries the fidelity function for each 
variation in the control values comes from the observa- 
tion that the forward evolution U {t, 0) and the backwards 
evolution U{T, t) need only be calculated once for a given 
control configuration. This in turn allows each derivative 
to be calculated with a constant number of matrix mul- 
tiplications, as opposed to the N required for the full 
Schrodinger evolution, leading to an overall scaling of 
0{N) as opposed to 0{N'^). 

The exact form of the derivative of the unitary time 
slice can be found in the original GRAPE paper, 



= -iAtHk.jUj, 

(3-7) 



1 

' " At 



For fine-grained control fields, i.e. ||-ffAt|| <C 1, the 
derivative can be approximated as 



duk,j 



-iAtHkUj. 



(3.8) 



This approximation can lead to difficulties in finding opti- 
mal solutions if the time step At is not sufficiently small, 
but this problem can be circumvented by evaluating the 
integral in Eq. (3.7) exactly [16]. Calculating this integral 
using the given form of Uj is straightforward after obtain- 
ing the diagonalization of H{jAt) (which is the preferred 
method of exponentiating the Hamiltonian given its her- 
miticity [35]). Thus, the form of the derivative of the time 
slice Uj with respect to the A;-th control [as in Eq. (3.7)] 
in the diagonal basis of the full Hamiltonian is 



OU ■ 6 ^ — 6 ^ 

K-|^|n,) = (m,|if.|n,)^— — ^ (3-9) 



for rij ^ rrij, or 

dU 

(mj|- — = —iAt{mj\Hk\'mj)e~ 

OUkJ 



(3.10) 



if Uj = m,j. Here, Am and A„. are eigenvalues of the full 
Hamiltonian at time slice j, associated with eigenvectors 
|TOj) and \nj) respectively. The actual derivative is then 
calculated by reverting the matrix elements back to the 
original, non-diagonal frame. 

The gradient formula, Eq. (3.5), explicitly assumes 
that the total Hamiltonian, control plus drift, is constant 
across each slice. For situations where the pulse shaping 
resolution limited by the pulse generator's sampling rate 
is slower than the Hamiltonian dynamics, it is important 
to consolidate the two time-scales. Several approaches 
have been tried that avoid changing the existing frame- 
work. The first is to discretize the dynamics into smaller 
slices and to penalize the differences between adjacent 
control values. This would lead to a search which has a 
dimension that scales with the smallest time-scale. An- 
other approach would be to coarse-grain in a way that 
completely discards the fast dynamics. For example, one 
can ignore the counter-rotating terms in the rotating 
wave approximation (RWA). The main purpose of the 
rest of the paper is to develop a method that redresses 
this tradeoff. 



IV. THE NEW ALGORITHM 

Our method is to separate the control parameters from 
the integration steps. This, for example, occurs when 
the control field arc combinations of Fourier components 
[17], when a control field arising from an arbitrary wave- 
form generator (AWG) is mixed with a carrier field of 
frequency much larger than the allowed sample rate of 
the AWG or when the control parameters are smoothed 
by a filter. In general, the field Ck{t) can depend on all 
input controls Ufc, though in many practical situations 
it only depends on local values of u^. Furthermore, our 
drift and control Hamiltonians can be time-dependent. 

As an example, we have conducted a "theorist's exper- 
iment" to show the typical response function given by an 
AWG to a set of digital inputs. In Fig. 1, the dashed 
purple line was inputed into a Tektronix AWG5014 (the 
AWG) and the response (solid green line) was measured 
on a LeCroy WavePro 725 Zi oscilloscope. This oscil- 
loscope had a bandwidth of 2.5 GHz (-3dB) which was 
large enough to ensure that all smoothing of the pulse 
was from the AWG with a 250 MHz bandwidth (-3dB). 
The pixels were set to Ins and 2ns in Fig. 1 (a) and 
(b), respectively. Here it is abundantly clear that the in- 
put control field and the waveform produced by current 
state-of-the-art AWG are not the same. The remaining 
two lines (red dot-dashed and blue dotted line) represent 
approximations that will be discussed in Sec. V. 

To numerically integrate the Schrodinger equation, it 
is a necessity to discretize the continuous fields Ck{t) into 
Sk,i, which is independented of our choice of u^. That 
is, time is digitized to t = I5t, where St is the time scale 
chosen such that ||^|| ||f|| < 1 and M6t = T. 
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FIG. 1: (color online) The purple dashed line represents the 
input controls Uj^k for a arbitrary chosen 5 pixel pulse with 
(a) Ins and (b) 2ns pixels, the green solid line is the output 
from the Tektronix AWG5014 running at 1 GSample/s (Ins 
pixels) and 500 MSamples/s (2 ns pixels) respectively with 
14 bit vertical resolution, the red dot-dashed line is the cubic 
spline interpolation, and the blue dotted line is the Gaussian 
filter approximation (see text for details). 



The evolution takes the form of Eq. (3.3) where now 





with the propagator. 



Ui = exp 



St 



(4.1) 



(4.2) 



Here, Hq j and j are time-sliced versions of HQ(t) and 
Hk{t). The update rule for uj, is computed the same as 
Eq. (2.4), but now the gradient of $ is found through 



duk. 



A/-1 

E 

1=0 



duk.-, dski ' 



(4.3) 



with 



a$ _ 2 



Re 



Tr U; 



rt 



ideal 



dU{T) 

'~a ' 

OSk,i 



Tr (C/idcalPQC/(r)tPQ) 



(4.4) 



where 

dU{T) 
dsk,i 



i+i 



n u„ 



=M-1 



dUi 
ds 



k,l 



n Un 



(4.5) 



Calculating the dUi/dsk^i proceeds exactly the same as 
in Sec. Ill, either exactly or through a linear approxi- 
mation. While this method will only ever be a linear 
approximation to the physical control fields, using the 
exact derivative as opposed to the linear approximation 
may speed up the gradient search. The partial derivative 
dsk,i/duk,j is provided by the transfer function Sk,i{uk)- 



V. TRANSFER FUNCTIONS 

There are a wealth of transfer functions for taking the 
controls U/^ to the continuous fields Ck{t). These transfer 
functions do not have to be linear but in many circum- 
stances a linear approximation is valid. For example, the 
filtering seen in Fig. 1 can is well approximated by a 
linear transfer function. 

A linear transfer function can be represented by 



JV-l 



Sk,l 



J' 



(5.1) 



3=0 



where the entirety of our transfer function is encapsu- 
lated in a transfer matrix, Tk^i.j. In this case, evaluating 
Eq. 4.3 is particularly straightforward since 



duk^j 



k,l,j- 



(5.2) 



In the case where the input controls are indexed 
by time, the transfer function (such as an interpolat- 
ing function, a filter, or a carrier function) will simply 
reparametrize the weights of the individual control pixels. 
If the transfer function provides variations on a smaller 
time scale than the original input pixels then additional 
sub-pixels Sk.i are required for the digitization of the con- 
tinuous field to be valid. This idea is illustrated in Fig. 
2, where in particular we can think of the gradient of 
the (purple-dashed) control pixels as the weighted sum of 
the gradients (the black arrows) of the sub-pixels (orange 
bars) inside or near the control pixel but which we cannot 
directly control. Note that this visualization works only 
when the control pixels are indexed in time and would 
not work e.g. if the controls were Fourier components. 



A. Cubic spline interpolation 

Lacking any detailed description of the transfer, 
smoothing can be approximated by interpolating between 
a discrete set of points with a piecewise cubic spline in- 
terpolation. Like the piecewise approximation in Eq. 3.1, 
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The non-zero matrix elements of Tfe,; j are then 
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FIG. 2: (color online) The purple dashed line represents the 
input controls pixels Uj,fc, the solid blue line represents the 
continuous field Cfc(t), and the orange bars represent the sub- 
pixel approximation (sk.j.i) to the continuous field. The black 
arrows represent the update vector of the individual sub-pixels 
while the large purple arrows represent the update (weighted 
average) for the control pixels. 



here we approximate the continuous fields as a piecewise 
continuous function 



N-l 

Ck{t) = ^f^'i^*^ (5-3) 

but now instead of remaining constant over the At inter- 
val the field is described by a cubic function, 



(5.4) 



To ensure that the field and its first derivative are contin- 
uous and to introduce the control parameters, we enforce 

Sk,j{{j + l/2)At) =uk,j, 



Uk.j + 1 — Uk^j-l 

' 2At '■ 



Sk,jiij + 3/2)At) =Uk,j+i, 

Uk.j+2 — UkJ 



(5.5) 



2At 



The box functions are offset by At/ 2 so that the con- 
trol parameters indicate the value at the center of the 
step and we additionally require that the function and 
its derivatives are zero at the boundaries (which can be 
enforced by padding the control vector with zeros). 

We can now derive the transfer matrix Tk^ij from the 
above conditions, and it is sparse. For each fc, I pair there 
are only four non-zero transfer matrix elements. There 
exists a / such that (/ -f 1/2) At < ISt < {f + 3/2) At. 



k,l,j' 



Tk,i,j' = 1 

T 

" 2At 



( 1)2 



ill 
2At3 

ill 
2At2 

^3 



ill 

2At2 

ili 

2At3 

^2 



(5.6) 



Ti 



k,l,j'+2 



2At3 2At2 



where t = ISt — [j' + 1/2) At. Spline fits of this form 
are shown in Fig. 1 as the red dot-dashed line. In these 
figures the spline is a better approximation for the 2ns 
pulse(Fig. lb) than for the Ins pulse (Fig. la). 



B. General Filters 

Most electronic systems undergo some amount of fil- 
tering. It is important to model this to understand the 
shape the waveforms will take upon reaching our quan- 
tum systems. To represent the effect of the filter the 
continuous field is 



/oo 
fk{t-t')uk{t')dt' 
-OO 

Cfe(t) =F-^[Fk{uj)F[ukm. 



(5.7) 



where fk{t), and Fk{uj) are the response functions of the 
filter for control k in the time and frequency domains 
respectively and J- is the Fourier transform. To ensure 
that the pulse is zero at the boundaries we add a finite 
number of control parameters to the start and end of the 
pulse that are fixed at zero. The number of parameters 
necessary for such a padding, n^, depends on the band- 
width of the filter, ujb, according to rir = \2Tr/ (w^ At)] . 
In the frequency domain, the transfer function is 



Sk,l{t) = ^ / Ukit') 



Ffc(w)e^"('**-*')dcjdt' (5.8) 



Using this with the piecewise constant control fields in 
Eq. (3.1) and the assumption that the filter function in 
the frequency domain is even, the transfer matrix be- 
comes 



'^kd,j - 



cos[w(Mt - (j + i)At)] sin[iwAt] 
Fk (uj) doj . 

, TTU! 

(5.9) 

It is easy to precompute Tkjj numerically, and for many 
filters it can be calculated analytically. Note also if — 
j| > rir we can assume that Tk.ij ~ 0. 

We demonstrate a calculation of Tk^ij using a Gaussian 
filter. In this case, the filter function is given by 



Fk{uj) = exp(-wVt^oJ 



(5.10) 
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where ojQ^ is the reference frequency [36] for the fc-th 
control. Performing the integration we find 



1 



erf 



erf 



Wo, 



^0. 



V 2 

ISt - {] + l)At 



(5.11) 



Gaussian filters are often used to approximate the ac- 
tual hardware filtering typically found in experiments, 
which are usually parametrized by their bandwidth ub 
(frequency of 3dB attenuation). For a Gaussian the 3dB 
attenuation occurs at lob = O.SSSTwq. Thus for the 
AWG with a bandwidth of u;B/27r = 250 MHz we find 
oJo/2-K = 425.4 MHz. The Gaussian approximation to the 
Tektronix AWG5014 is shown in Fig. 1 as the blue dotted 
line. Here, we see that it reasonably well approximates 
the experiment (much better then the spline) and for the 
remainder of the paper, we used this as our benchmark. 



VI. APPLICATIONS 

To demonstrate the usefulness of this approach, we 
consider three examples in this sections. The first is the 
implementation of a tt pulse in an anharmonic oscillator 
where the RWA approximation has been made. This sys- 
tem has been studied in detail in Ref. [10, 26, 32, 38, 39] 
and is very applicable to the superconducting qubit com- 
munity. Here, we enforce that the controls have limited 
bandwidth (we choose a bandwidth consistent with state- 
of-the-art electronics) and find optimal pulses with our 
filtering and spline smoothing techniques. The second ex- 
ample is the same system but without the RWA. This al- 
lows us to demonstrate the carrier modulation technique 
as well as robustness of the fidelity to the initial rela- 
tive phase between the carrier and the shaped envelope. 
The last example is a \/iSWAP between two off-resonant 
qubits coupled by a XX + YY interaction. This is a 
demonstration of multiple tones and a time-dependent 
drift used to couple the off-resonant qubits. 



C. Carrier modulation 

In the above, we consider transfer functions that 
smooth a piecewise function into a continuous function. 
However, there is an interesting situation that arises 
when the control fields are mixed with a carrier. In this 
case, the transfer matrix is 

Tk,i,j = f{l5t)nj{m,At), (5.12) 

where f{t) is the carrier function which in most cases 
takes the form of a cosine function with an arbitrary un- 
known reference phase tp. In appendix A and Sec. VI B 
we show how to adjust the cost function $ to be robust 
against this uncertainty. Note we could also easily com- 
bine the previous smoothing transfer functions with this 
carrier effect. 



D. Fourier components 

In previous sections, we envision the control parame- 
ters as the magnitudes of piecewise constant functions, 
i.e. square pulses, and the approximations as smoothed 
versions of these fields. In principle, we could instead 
look at the Wfe^'s as more arbitrary parameters such as 
Fourier components (see e.g. Ref. [37]), where now we 
write 

Cfe(*) = X] sin(wjt), (5.13) 

3 

and the transfer matrix is simply 

Tk,i,j ~ sm.{(jjjl6t). (5.14) 



A. Optimization with large filtering effects 

For the first example, we consider a system with the 
following Hamiltonian: 

H = [S''{t) cos{ujit + V) + £^it) sm{uit + V)] (F + F^) 

(6.1) 

where ioi is the — 1 transition frequency, 1^2 is the — 
2 transition frequency, ip is the unknown phase of the 
carrier signal at the start of the pulse, F is the effective 
lowering operator F = |0)(1| -|- ■\/2|l)(2|, and £x{t) and 
Sy{t) are the controls. This system represents the first 
three levels of an anharmonic oscillator and effectively 
models many systems for simulating a qubit. With this 
system we want to perform a tt rotation in the subspace 
formed by the and 1 states without losing population to 
the third level. It is standard to define this n pulse in a 
frame rotating at the frequency uj\ . Moving to this frame 
and making the rotating wave approximation (RWA), we 
can model this system by the Hamiltonian 

H = g^(^)£±l!l + sy(t) ~ '^^^ + A|2)(2|, (6.2) 

where A = 2wi — u!2 ■ 

Under the RWA with two controls (with infinite band- 
width and sampling rate) it is always possible to find 
a solution that gives a perfect tt pulse [10]. However, 
when we restrict the control parameters to Ins (a typical 
setting for current microwave AWGs), we find that we 
need a gate time of least 4ns (4 pixels) to reach errors 
below 10~^. This is shown in Fig. 3 as the dotted purple 
line where we have taken A/27r = —500MHz . Naively, 
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FIG. 3: (colour online) Gate error as a function of the gate 
time for the filtering example. The dotted purple line is the 
optimal grape solution without any filtering. The dashed 
green line is the predicted error when a the optimal grape 
solutions are filtered by a Gaussian filter approximating a 
AWG with a 250 MHz bandwidth. The dot-dashed cyan line 
is the predicted gate error after filtering with the Gaussian 
filter for the spline optimization with 2 sub-pixels per control 
pixel. The solid red line and the dotted blue line are the pre- 
dicted gate errors after filtering with the Gaussian filter for 
the Gaussian filter optimization with 2 and 20 sub-pixels per 
control pixel, respectively. The vertical dashed line indicates 
the gate time used in Fig. 4. Other parameters are given in 
the text. 



one would then predict that these optimized pulses could 
be perfectly implemented. However, due to the internal 
filtering imposed by the AWG this is not the case (see 
Fig. 1). To demonstrate this we take the numerically 
optimized pulses and filter them with a Gaussian filter 
approximation to a AWG with an internal bandwidth of 
250MHz (see Sec. VB), which for short we will call the 
'AWG filter'. The predicted error is then shown in Fig 
3 as the dotted green line. Here we see that the effect 
of the filtering is drastic, and hence the optimized pulses 
will not perform well for quantum operations with the 
AWG. 

Taking the filtering into account during the optimiza- 
tion allows for much better pulses. We can do this either 
by finding smooth pulses with a cubic spline interpola- 
tion or by actually taking into account the filter. To 
demonstrate this we plot in Fig. 3 the error as a func- 
tion of gate time when optimizing under the cubic spline 
interpolation (dot-dashed cyan line) and a Gaussian filter 
(solid red line). Here we included 2 subpixels per control 
pixel, and after finding the optimized solutions, applied 
the AWG filter. Increasing the number of subdivisions 
allows our algorithm to better approximate the AWG fil- 
ter; this is illustrated by the dotted blue line where we 
find for 20 sub-pixels a greatly improved performance for 
all gate times. To get an indication of the performance 
of our algorithm with the number of subdivisions, we set 
the gate time to 4ns (vertical dotted line in Fig. 3) and 
plot in Fig. 4 the gate error as a function of number of 
sub-pixels. Here, we observed that for only a few subdivi- 
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FIG. 4: (colour online) Gate error as a function of the number 
of subpixels for the filtering example. The dotted purple line 
is the optimal grape solution without any filtering (limited 
by the numerical precission). The dashed green line is the 
filtered optimal grape solutions when the filter is a Gaussian 
filter approximating a AWG with a 250 MHz bandwidth. The 
dot-dashed cyan and solid red lines are the predicted gate 
errors after filtering with the Gaussian filter for the spline and 
Gaussian filter optimizations, respectively. Other parameters 
are given in text. 

sions the performance of our algorithm reaches very small 
error rates. We also find that the spline optimization is 
not as good as the Gaussian filter. This is expected as 
we have assumed the real situation (AWG filter) to be 
a Gaussian filter, hence the spline optimization is not 
necessarily going to find pulses that are consistent with 
the AWG filter. On the other hand, picking the correct 
transfer function improves the situation, even for very 
few subdivisions. 



B. Optimization with time-dependent carriers 

For the second example, we take the same system and 
goal as the first example, but we optimize without mak- 
ing the RWA. Instead of focusing on the errors the arise 
from smoothing, we are concerned with the errors aris- 
ing from the rotating terms in order to demonstrate the 
carrier modulation optimization. In this case, the Hamil- 
tonian is 

\T(1 + p-Sif^it-ai-iA) + h cl 

( C\ 

, [zr(l + e-2^"i*-2^^)+h.c] 

2 ■ 

To demonstrate the error from using the wrong Hamil- 
tonian, we first find optimal solutions for the tt pulses 
assuming the RWA [optimization with Eq. (6.2)] for 
A/27r = -500 MHz, a;i/27r = 2.0 GHz, and control pix- 
els of 0.125 ns. This optimal solution is shown in Fig. 
5 as the solid purple line as a function of the gate time. 
For each gate time the controls are then used to evalu- 
ate what the fidelity of the operation would be if we did 
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FIG. 5: (colour online) Gate error as a function of gate time 
for the carrier example. The horizontal solid purple line 
is the (target and achieved) gate error for the optimal so- 
lution when we make the RWA. The dotted blue lines are 
the predicted error when we do not make the RWA but find 
the optimal solution with the RWA. The lower dashed green 
lines are the gate errors when we optimize without making 
the RWA. Each line represents one of the nine different val- 
ues of the relative phase between the envelope and the car- 
rier: 2.46245, 2.13875, 1.57081, 0.304685, 0.043838, 1.65238, 
0.914728, 2.02047, 0.518253. 



not make the RWA for various phases V' [evolution un- 
der Eq. (6.3)]. These results are shown in Fig. 5 as the 
dotted blue lines, where each line represents a randomly 
chosen phase. This figure clearly indicates that for short 
gate times, neglecting the rotating terms leads to a large 
error. 

To perform the optimization with the rotating terms 
we simply include them as a carrier function. Since we 
would also like to find pulses which don't depend on ip, 
we also use the robust technique outlined in Appendix 
A. This is essentially an optimization of average fidelity 
over all possible ip. In Fig. 5 we plot the gate error 
as a function of gate time (dashed green lines), seeing 
that it is possible to find pulses that robustly remove the 
rotation errors. Here we set the number of sub-pixels per 
control pixel to 100. 



C. Optimization with multiple tones 



FIG. 6: (colour online) Gate error as a function of the number 
of subpixels for the multi-tone example. The solid red line is 
for A/2-K = 1.0 GHz and the green dashed line is for A/2tt = 
0.5 GHz. Other parameters are given in the text. 



the form, 

H =£P (t)4^) cos (c. Wt) + £W 

+ (t)a(^) cos (J^k) + £f ) (i)ai^) sin (.ft) 
+ J + ai^Vf ) + u^.a^^^ + u;,ai'\ 

(6.4) 

where the superscripts index the qubit, a± are the raising 
and lowering operators, J is the strength of the coupling, 
uji and UJ2 are the qubit frequencies, and wf and wf 
are the drive frequencies. Choosing the frame rotating at 
the energies of the two qubits (which is conventional, but 
not the only available frame) and setting the first drive 
to wf = (^1+1^2 g^jj^^ -(.jjg second drive to wf = 012, we 
find 

+ J(aV')ai^)e^^*+ai^)ai^)e-^*), 

(6.5) 



The third example that we consider demonstrates that 
fast dynamics can occur simply from energy differences 
between different components in a system. We consider 
two off-resonant qubits coupled by a XX + YY inter- 
action and drive with two different frequencies, one res- 
onant with one of the qubit 's transition frequency and 
the second tuned to exactly the average of the qubits' 
transition frequencies, in order to generate an VlSWAP 
operation. Thus, we see the system has three frequen- 
cies, and there is no frame in which all fast dynamics can 
disappear. Specifically, we start with the Hamiltonian of 



where A = wi — 1^2 is the energy difference between the 
qubits. 

To illustrate the scaling with subpixels for this exam- 
ple we choose a gate time of tg = 20ns, J/27r — 94 MHz, 
and A/27r = 0.5 or 1.0, and in Fig. 6 we plot the pre- 
dicted error as a function of the number of sub-pixels 
for Ins pixels. This error is calculated by integrating 
the Schrodinger equation for the optimized controls on a 
much finer grid. As expected, we see that for the larger 
detuning much more subpixels are required to achieve the 
same low error. 
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VII. CONCLUSIONS 

In this article, we have developed a numerical method 
for fincling optimal control pulses for implementing uni- 
tary gates when the control variation is at a time scale 
that is much different than the numerical integration time 
scale (set by the time variation of the Hamiltonian) . This 
occurs in many physical examples ranging from the case 
when the controls are variations of an envelope function 
on a carrier field to the case when the controls are fil- 
tered. We have shown that it is necessary to separate 
the input controls (occurring at the sampling rate of the 
waveform generator) from the fast-varying dynamics of 
the Hamiltonian. We accomplish this by distinguishing 
the control discretization from the discretization neces- 
sary for integration. We have shown that these dynamics 
can be computed quickly. Furthermore, we have given 
some examples of how these techniques can be applied to 
relevant physical systems. 
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Appendix A: Robust controls 

In many cases, we want to find controls that are robust 
to some uncertain parameter [40-42]. This could be, 
for example, an uncontrolled phase between a carrier and 
the the control pixels. To account for this, we define an 
average performance index 

$ = / (Al) 

where is defined the same way as Eq. (2.3) for each 
possible tp. Since this is a linear combination of all the 
possible the gradient of $ is found by 

= / ^di^- A2 

dukj duk,j 

In practice the integral is replaced by summation, and 
we sample discretely over all possible tp's. 
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